1 Import libraries

source("R/bioage_estimate_median.R")
source("R/gompertz_draw.R")
source("R/weibull_draw.R")
source("R/generate_data.R")
library(tidyverse)
library(MASS)
library(psbcGroup)
library(brms)
library(glue)
library(bayesplot)

library(cmdstanr) # Faster STAN

2 Generation of survival data

generate_data = function(pop_n, obs_n, p, g, ltname, seed, force_recalc, X_plots, 
                         X_rho, β_rho, non_zero_betas, X_scale, active_hazard, 
                         betasep, bscale, betaplot = F, target_snr = 2, 
                         method = "weibull", a, b) {
  
  print("Generating βs...")
  betas = generate_betas(p = p, g=g, rho = β_rho, rho_between = 0,
                       seed = seed, mu_u = betasep, mu_l = -betasep, 
                       beta_scale = bscale, plot = betaplot, 
                       non_zero_groups = non_zero_betas, 
                       active_hazard = active_hazard,
                       target_snr = target_snr,
                       method = method, # Different scaling for different methods
                       weib_shape = a, weib_scale = b,
                       gomp_a = a, gomp_b = b)
  
  # betas$beta = betas$beta + 0.01
  
  
  print("Generating X...")
  X = generate_X(n = pop_n, p = p, g = g, rho = X_rho, seed = seed,
                 rho_between = 0, X_plots = X_plots, scale = X_scale)
  print("Done!")

  
  if (method == "weibull") {
    # a = shape
    # b = scale
    print("Generating Weibull population lifetable...")
    lt = generate_population_lifetable_weibull(N_pop = pop_n, M = p, 
                                              betas = betas$beta, X = X$X, 
                                              filename = ltname, seed = seed, 
                                              force_recalc = force_recalc, 
                                              a = a, b = b)  # Weibull params
    print("Done!")
    
    plot = ggplot(lt, aes(x = t, y = mrl)) +
      geom_line() +
      labs(title = "Weibull Population Lifetable MRL",
           x = "Chronological age (years)", y = "Mean Residual Life") +
      theme_minimal()
    print(plot)
    
    print("Generating data from Weibull lifetable...")
    df_sim = create_dataset_weibull(M = p, n_obs = obs_n, G = g,
                                   gsize = (p/g), lt = lt, betas = betas$beta,
                                   seednr = seed+1, X_plots = F, 
                                   a = a, b = b,  # Weibull params
                                   X_rho = X_rho, X_scale = X_scale, 
                                   followup = 20) # Nr years followed
    
  } else {
    print("Generating Gompertz population lifetable...")
    lt = generate_population_lifetable_gompertz(N_pop = pop_n, M = p, 
                                               betas = betas$beta, X = X$X, 
                                               filename = ltname, seed = seed, 
                                               force_recalc = force_recalc, 
                                               a = a, b = b)  # Gompertz params
    print("Done!")
    
    plot = ggplot(lt, aes(x = t, y = mrl)) +
      geom_line() +
      labs(title = "Gompertz Population Lifetable MRL",
           x = "Chronological age (years)", y = "Mean Residual Life") +
      theme_minimal()
    print(plot)
    
    print("Generating data from Gompertz lifetable...")
    df_sim = create_dataset_gompretz(M = p, n_obs = obs_n, G = g,
                                    gsize = (p/g), lt = lt, betas = betas$beta,
                                    seednr = seed+1, X_plots = F, 
                                    a = a, b = b,  # Gompertz params
                                    X_rho = X_rho, X_scale = X_scale, 
                                    followup = 20) # Nr years followed
  }
  
  print("Done!")
  return(list(df = df_sim, true_betas = betas, lt = lt))
}

3 Functions to estimate MRL

3.1 Frequentist AFT

aft_reg = function(df_sim, p, method) {
  result = switch(method,
    "gompertz" = aft_reg_gompertz(df_sim, p),
    "weibull" = aft_reg_weibull(df_sim, p),
    stop("Unsupported method: ", method)
  )
  return(result)
}

3.1.1 Gompertz implementation

aft_reg_gompertz = function(df_sim, p) {
  true_betas = df_sim$true_betas
  lt = df_sim$lt
  df_sim = df_sim$df

  
  # Train/test split 50/50
  train_indices = sample(1:nrow(df_sim), nrow(df_sim)/2)
  df_sim_train = df_sim[train_indices,]
  df_sim_test = df_sim[-train_indices,]
  
  
  # Does not work in highdim cases
  if (dim(df_sim_train)[1] < p) {
    return(NULL)
  }
  
  
  #### FIT MODEL
  predictor_vars = paste(c(glue("x{1:p}")), collapse = " + ")
  aft_formula = as.formula(paste("Surv(age_start, age_end, status) ~", predictor_vars))
  
  fit_aft = tryCatch({
    eha::aftreg(
      formula = aft_formula, 
      data = df_sim_train, 
      dist = "gompertz",
      param = "lifeExp",
      control = list(maxit = 100, trace = FALSE)
    )
  }, error = function(e) {
    message("AFT model error: ", e$message)
    NULL
  })
  
  if (is.null(fit_aft)) {
    print("Fit failed")
    return(NULL)
  }

  sigma_est <- exp(fit_aft$coefficients["log(scale)"])
  tau_est <- exp(fit_aft$coefficients["log(shape)"])
  a_est <- tau_est / sigma_est
  b_est <- 1 / sigma_est
  
  betas_est_aft <- fit_aft$coefficients[1:p]
  linpred_aft <- rowSums(sweep(df_sim_test[,1:p], 2, betas_est_aft, "*"))
  
  distr_ests = list(sigma = sigma_est, tau = tau_est, a = a_est, b = b_est)
  
  #### ESTIMATE MRL
  df_sim_test$pred_mrl = NA
  
  for (i in 1:nrow(df_sim_test)) {
    mrl_uncon <- integrate(gomp_baseline_surv, lower = (df_sim_test$age_start[i] * exp(linpred_aft[i])), 
                         upper=Inf, a = a_est, b = b_est)$value
    s_cond <-  gomp_baseline_surv(df_sim_test$age_start[i] * exp(linpred_aft[i]), a = a_est, b = b_est) 
    df_sim_test$pred_mrl[i] <- (mrl_uncon / s_cond) * exp(-linpred_aft[i])
    
    if(is.na(df_sim_test$pred_mrl[i]) || df_sim_test$pred_mrl[i] < 0) {
      df_sim_test$pred_mrl[i] <- NA
    }
  }
  
  pred_mrl = df_sim_test$pred_mrl
  # obtain biological age, assuming true lifetable is known
  aft_b = vector(length = length(pred_mrl))
  
  for (i in 1:length(pred_mrl)){
    aft_b[i] = lt$t[which.min(abs(lt$mrl - pred_mrl[i]))]
  }

  # Calculate rmses
  aftreg_coef_rmse = sqrt(mean((betas_est_aft - true_betas$beta)^2))
  aftreg_rmse = sqrt(mean((df_sim_test$mrl - df_sim_test$pred_mrl)^2))
  aftreg_bioage_rmse = sqrt(mean((aft_b - df_sim_test$b_age)^2))

  # plt = df_sim_test %>% ggplot(aes(x = mrl, y = pred_mrl)) +
  #   geom_point() +
  #   geom_abline(color = 'red') +
  #   labs(title = "True vs Predicted MRL for EHA AFT Gompertz",
  #        subtitle = glue("RMSE = {aftreg_rmse}")) +
  #   theme_minimal()
  # print(plt)

  true_betas$betahat = betas_est_aft

  # plt = true_betas %>%
  #   mutate(index = 1:p) %>%
  #   ggplot(aes(x = index)) +
  #   geom_segment(aes(xend = index, y = beta, yend = betahat), color = "gray") +
  #   geom_point(aes(y = beta, color = "True")) +
  #   geom_point(aes(y = betahat, color = "Predicted")) +
  #   scale_color_manual(values = c("True" = "black", "Predicted" = "skyblue"),
  #                      name = "Type") +
  #   ylab("β") +
  #   xlab("Index") +
  #   labs(title = "True vs Predicted β-coefficients for AFT Gompertz", 
  #        subtitle = glue("RMSE = {aftreg_coef_rmse %>% round(4)}")) +
  #   theme_minimal()
  # print(plt)
  
  return(list(
    test_data = df_sim_test,
    true_betas = true_betas,
    predicted_betas = betas_est_aft,
    predicted_bioage = aft_b,
    rmses = list(coef = aftreg_coef_rmse, mrl = aftreg_rmse, bioage = aftreg_bioage_rmse),
    method = "AFT gompertz",
    distr_ests = distr_ests
  ))
}

3.1.2 Weibull implementation

aft_reg_weibull = function(df_sim, p) {
  true_betas = df_sim$true_betas
  lt = df_sim$lt
  df_sim = df_sim$df

  
  # Train/test split 50/50
  train_indices = sample(1:nrow(df_sim), nrow(df_sim)/2)
  df_sim_train = df_sim[train_indices,]
  df_sim_test = df_sim[-train_indices,]
  
  # Does not work in highdim cases
  if (dim(df_sim_train)[1] < p) {
    return(NULL)
  }
  
  #### FIT MODEL
  predictor_vars = paste(c(glue("x{1:p}")), collapse = " + ")
  aft_formula = as.formula(paste("Surv(age_start, age_end, status) ~", predictor_vars))
  
  fit_aft = tryCatch({
    eha::aftreg(
      formula = aft_formula, 
      data = df_sim_train, 
      dist = "weibull",
      param = "lifeExp", # Otherwise beta signs are inverted
      control = list(maxit = 100, trace = FALSE)
    )
  }, error = function(e) {
    message("AFT model error: ", e$message)
    NULL
  })
  fa <<- fit_aft
  if (is.null(fit_aft)) {
    print("Fit failed")
    return(NULL)
  }
  print("----------")
  shape_est <- exp(fit_aft$coefficients["log(shape)"])
  cat("shape_est:", shape_est, "\n")
  scale_est <- exp(fit_aft$coefficients["log(scale)"])
  cat("scale_est:", scale_est, "\n")
  
  betas_est_aft <- fit_aft$coefficients[1:p]
  linpred_aft <- rowSums(sweep(df_sim_test[,1:p], 2, betas_est_aft, "*"))
  
  
  #  scale^(-shape) == lambda so we use it as est (from Marije's script)
  lambda_est <- scale_est^(-shape_est)
  nu_est <- shape_est
  
  #### ESTIMATE MRL
  df_sim_test$pred_mrl = NA
  
  for (i in 1:nrow(df_sim_test)) {
    scale_adj = scale_est * exp(linpred_aft[i])

    mrl_uncon <- integrate(function(t) {
      1 - pweibull(t, shape = shape_est, scale = scale_adj)
    }, lower = df_sim_test$age_start[i], upper = Inf)$value

    s_cond <- 1 - pweibull(df_sim_test$age_start[i], shape = shape_est, scale = scale_adj)
    
    df_sim_test$pred_mrl[i] <- mrl_uncon / s_cond

    # mrl_uncon <- integrate(weib_baseline_surv, 
    #                       lower = (df_sim_test$age_start[i] * exp(linpred_aft[i])), 
    #                       upper = Inf, 
    #                       lambda = lambda_est, nu = nu_est)$value
    # s_cond <-  weib_baseline_surv(df_sim_test$age_start[i] * exp(linpred_aft[i]), 
    #                              lambda = lambda_est, nu = nu_est) 
    # df_sim_test$pred_mrl[i] <- (mrl_uncon / s_cond) * exp(-linpred_aft[i])
  }

  print("----------")
  # obtain biological age, assuming true lifetable is known
  pred_mrl = df_sim_test$pred_mrl
  aft_b = vector(length = length(pred_mrl))
  
  for (i in 1:length(pred_mrl)){
    aft_b[i] = lt$t[which.min(abs(lt$mrl - pred_mrl[i]))]
  }
  
  # Calculate rmses
  aftreg_coef_rmse = sqrt(mean((betas_est_aft - true_betas$beta)^2))
  aftreg_rmse = sqrt(mean((df_sim_test$mrl - df_sim_test$pred_mrl)^2))
  aftreg_bioage_rmse = sqrt(mean((aft_b - df_sim_test$b_age)^2))
  
  plt = df_sim_test %>% ggplot(aes(x = mrl, y = pred_mrl)) +
    geom_point() +
    geom_abline(color = 'red') +
    labs(title = "True vs Predicted MRL for EHA AFT Weibull",
         subtitle = glue("RMSE = {aftreg_rmse}")) +
    theme_minimal()
  # print(plt)
  
  true_betas$betahat = betas_est_aft
  
  plt = true_betas %>%
    mutate(index = 1:p) %>%
    ggplot(aes(x = index)) +
    geom_segment(aes(xend = index, y = beta, yend = betahat), color = "gray") +
    geom_point(aes(y = beta, color = "True")) +
    geom_point(aes(y = betahat, color = "Predicted")) +
    scale_color_manual(values = c("True" = "black", "Predicted" = "skyblue"),
                       name = "Type") +
    ylab("β") +
    xlab("Index") +
    labs(title = "True vs Predicted β-coefficients for AFT Weibull", 
         subtitle = glue("RMSE = {aftreg_coef_rmse %>% round(4)}")) +
    theme_minimal()
  # print(plt)
  
  return(list(
    test_data = df_sim_test,
    true_betas = true_betas,
    predicted_betas = betas_est_aft,
    predicted_bioage = aft_b,
    rmses = list(coef = aftreg_coef_rmse, mrl = aftreg_rmse, bioage = aftreg_bioage_rmse),
    method = "AFT weibull",
    distr_ests = list(shape = shape_est, scale = scale_est)
  ))
}

3.2 PSBC Bayesian AFT

psbc_reg = function(df_sim, p, g) {
  true_betas = df_sim$true_betas
  df_sim = df_sim$df
  lt = df_sim$lt
  
  # Train/test split 50/50
  train_indices = sample(1:nrow(df_sim), nrow(df_sim)/2)
  
  df_sim_train = df_sim[train_indices,]
  df_sim_test = df_sim[-train_indices,]
  
  Y_train = cbind(df_sim_train$age_start, df_sim_train$age_end, df_sim_train$status)
  Y_test = cbind(df_sim_test$age_start, df_sim_test$age_end, df_sim_test$status)
  
  # Create covariate matrix X
  X_train = as.matrix(df_sim_train[, paste0("x", 1:p)])
  X_test = as.matrix(df_sim_test[, paste0("x", 1:p)])
  
  XC = NULL # CONFOUNDERS
  
  # Create group indicator for each variable
  grpInx = rep(1:g, each = p/g)
  
  # Set hyperparameters (PRIOR)?
  hyperParams = list(
    a.sigSq = 2, b.sigSq = 2,
    mu0 = 0, h0 = 10^6,
    v = 10^6
  )
  
  # Set starting values https://cran.r-project.org/web/packages/psbcGroup/psbcGroup.pdf
  startValues = list(
    w = log(Y_train[,1]),
    beta = runif(p, -0.5, 0.5), # uniform draws from realistic range
    tauSq = rep(0.1, p),
    mu = 0.1,
    sigSq = 0.1,
    lambdaSq = 0.1,
    betaC = rep(0.11, 0)
  )

  
  # MCMC parameters format __ CHANGE __ 
  mcmcParams = list(
    run = list(
      numReps = 5000,
      thin = 5, # Update params every n steps
      burninPerc = 0.2 # Discard the first ... steps
    ),
    tuning = list(
      mu.prop.var = 0.1,
      sigSq.prop.var = 0.005,
      L.beC = 50,
      M.beC = 1,
      eps.beC = 0.001,
      L.be = 100,
      M.be = 1,
      eps.be = 0.001
    )
  )
  
  fit_bayes_aft = aftGL_LT(Y_train, X_train, XC, grpInx, hyperParams, startValues, mcmcParams)
  beta_samples = fit_bayes_aft$beta.p

  # # Convert to long format
  # beta_long = beta_samples %>%
  #   as.data.frame() %>%
  #   mutate(iteration = row_number()) %>%
  #   pivot_longer(cols = -iteration, 
  #                names_to = "beta_index", 
  #                values_to = "value") %>%
  #   mutate(beta_index = factor(beta_index, levels = paste0("V", 1:ncol(beta_samples))))
  
  # # Create boxplot
  # plt = ggplot(beta_long, aes(x = beta_index, y = value)) +
  #   geom_boxplot(fill = "lightblue", alpha = 0.7) +
  #   geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  #   labs(
  #     title = "Posterior Distributions of Beta Coefficients",
  #     x = "Beta Index",
  #     y = "Coefficient Value"
  #   ) +
  #   theme_minimal() +
  # theme(axis.text.x = element_text(angle = 45, hjust = 1))
  # 
  # print(plt)

  
  # Extract betas using the maximum likelihood via density
  get_mode = function(x){
    xdist = density(x)
    mode = xdist$x[which.max(xdist$y)]
    return(mode)
  }
 # betas = fit_bayes_aft$beta.p %>% apply(MARGIN = 2, FUN = get_mode)
  betas = apply(fit_bayes_aft$beta.p, 2, mean)
  # Make linear predictor for MRL
  X_test = X_test %>% as.matrix
  linpred = X_test %*% betas
  
  # Change to lognormal
  # Create baseline lognormal survival
  # Implemented in R -> check log scales
  # dlnorm, use any of them 1 - F(x)
  # instead of gomp_baseline_surv no false, 1-
  
  
  # a = exp(-9)
  # b = 0.085
  
  # 1-plnorm(1:100, meanlog = fit_bayes_aft$mu.p[120], sdlog = fit_bayes_aft$sigSq.p[120])
  
  # Estimate mean from MCMC draws
  mu_est = mean(fit_bayes_aft$mu.p)
  sigSq_est = mean(fit_bayes_aft$sigSq.p)

  for (i in 1:nrow(df_sim_test)){
    mrl_uncon = integrate(function(x) 1 - plnorm(x, meanlog = mu_est, sdlog = sqrt(sigSq_est)), # Var to Sd via sqrt
                       lower = (df_sim_test$age_start[i] * exp(linpred[i])), 
                       upper = Inf)$value
    s_cond = 1 - plnorm(df_sim_test$age_start[i] * exp(linpred[i]), meanlog = mu_est, sdlog = sqrt(sigSq_est))
    df_sim_test$pred_mrl[i] <- (mrl_uncon / s_cond) * exp(-linpred[i])
  }
  
  # Get the RMSE
  pred_RMSE = sqrt(mean((df_sim_test$mrl - df_sim_test$pred_mrl)^2))

  beta_RMSE =sqrt(mean((betas - true_betas$beta)^2))

  # MRL prediction plot
  plt = df_sim_test %>% ggplot(aes(x = mrl, y = pred_mrl)) +
    geom_point() +
    geom_abline(color = 'red') +
    labs(title = "True vs Predicted MRL for PSBC Bayesian AFT",
         subtitle = glue("RMSE = {pred_RMSE}")) +
    theme_minimal()

  # print(plt)

  # Plot difference between true and estimates betas
  true_betas$betahat = betas

  plt = true_betas %>%
    mutate(index = 1:p) %>%
    ggplot(aes(x = index)) +
    # Vertical lines between pred and actual
    geom_segment(aes(xend = index, y = beta, yend = betahat), color = "gray") +
    geom_point(aes(y = beta, color = "True")) +
    geom_point(aes(y = betahat, color = "Predicted")) +
    scale_color_manual(values = c("True" = "black", "Predicted" = "skyblue"),
                       name = "Type") +
    ylab("β") +
    xlab("Index") +
    labs(title = "True vs Predicted β-coefficients for Bayesian AFT", subtitle = glue("RMSE = {beta_RMSE %>% round(4)}")) +
    theme_minimal()
  # print(plt)
  

  return(list(
    test_data = df_sim_test,
    true_betas = true_betas,
    predicted_betas = betas,
    method = "PSBC",
    distr_ests = list(mu = mu_est, sigSq = sigSq_est)
  ))
}

3.3 BRMS Bayesian AFT

brms_reg = function(df_sim, p, g, method, seed = 123) {
  result = switch(method,
    "gompertz" = brms_reg_gompertz(df_sim, p, g, seed),
    "weibull" = brms_reg_weibull(df_sim, p, g, seed),
    stop("Unsupported method: ", method)
  )
  return(result)
}

3.3.1 Gompertz implementation

brms_reg_gompertz = function(df_sim, p, g, seed) {
  # Grompertz stan implementation source:
  # https://metrumrg.com/wp-content/uploads/2024/12/2024-Yoder-cure-rate-brms.pdf
  # Todd Yoder, Andrew Tredennick, Timothy Waterhouse: Gompertz Cure Rate Survival Models with Stan and brms
  
  stan_funs <- "
    real gompertz_lpdf(real t) {
      return log(mu) + gamma*t + exp(gamma*t);
    }
    real gompertz_lcdf(real t) {
      return log(1 - exp((1 - exp(gamma*t))));
    }
    real gompertz_lccdf(real t) {
      return (1 - exp(gamma*t));
    }
    real gompertz_rng(real mu, real gamma) {
      real sim_time;
      sim_time = 1/gamma*log(1 - gamma/mu*log(1 - uniform_rng(0, 1)));
      if (is_nan(sim_time))
        sim_time = positive_infinity();
      return sim_time;
    }
  "
  
  gompertz_fam <- brms::custom_family(
    name = "gompertz",
    dpars = c("mu", "gamma"),
    links = c("log", "identity"),
    lb = c(0, NA),
    type = "real",
    log_lik = function(i, prep) {
      mu <- brms::get_dpar(prep, "mu", i = i)
      gamma <- brms::get_dpar(prep, "gamma", i = i)
      t <- prep$data$Y[i]
      cens <- prep$data$cens[i]
      if (cens == 0) x <- gompertz_lpdf(t, mu, gamma)
      if (cens == 1) x <- gompertz_lccdf(t, mu, gamma)
      return(x)
    },
    posterior_predict = function(i, prep, ...) {
      mu <- brms::get_dpar(prep, "mu", i = i)
      gamma <- brms::get_dpar(prep, "gamma", i = i)
      return(gompertz_rng(mu, gamma))
    }
  )
  
  true_betas = df_sim$true_betas
  lt = df_sim$lt
  df_sim = df_sim$df
  
  # Train/test split 50/50
  train_indices = sample(1:nrow(df_sim), nrow(df_sim)/2)
  df_sim_train = df_sim[train_indices,]
  df_sim_test = df_sim[-train_indices,]

  #### FIT MODEL
  preds = paste(paste0("x", 1:p), collapse = " + ")
  formula = as.formula(paste("age_end | cens(1 - status) + trunc(lb = age_start) ~", preds))

  priors = c(
    prior(normal(0, 0.5), class = "b"), # BETAS, using ridge, λ = 2 in this case
    prior(normal(log(80), 0.4), class = Intercept)#, # about human life expectancy
    # prior(gamma(2,1), class = shape) # Hazard
  )

  fit_bayes_aft = brm(
    formula = formula,
    data = df_sim_train,
    family = gompertz_fam,
    stanvars = stanvar(scode = stan_funs, block = "functions"),
    prior = priors,
    warmup = 2000, iter = 4000, chains = 4, cores = 4
  )
  
  # Extract posterior samples and calculate estimates
  post = posterior_samples(fit_bayes_aft)
  beta_cols = paste0("b_x", 1:p)
  betahat = apply(post[, beta_cols], 2, mean)
  
  X_test = df_sim_test[,1:p] %>% as.matrix
  linpred = X_test %*% betahat
  
  # Extract Gompertz parameters from posterior
  mu_est = mean(post$mu)
  gamma_est = mean(post$gamma)
  
  # Convert to a/b parameterization for consistency with AFT
  a_est = gamma_est
  b_est = 1/mu_est
  
  #### ESTIMATE MRL (using median approach like AFT)
  df_sim_test$pred_mrl = NA
  
  for (i in 1:nrow(df_sim_test)) {
    mrl_uncon <- integrate(gomp_baseline_surv, lower = (df_sim_test$age_start[i] * exp(linpred[i])), 
                         upper=Inf, a = a_est, b = b_est)$value
    s_cond <-  gomp_baseline_surv(df_sim_test$age_start[i] * exp(linpred[i]), a = a_est, b = b_est) 
    df_sim_test$pred_mrl[i] <- (mrl_uncon / s_cond) * exp(-linpred[i])
    
    if(is.na(df_sim_test$pred_mrl[i]) || df_sim_test$pred_mrl[i] < 0) {
      df_sim_test$pred_mrl[i] <- NA
    }
  }
  
  # obtain biological age, assuming true lifetable is known
  pred_mrl = df_sim_test$pred_mrl
  brms_b = vector(length = length(pred_mrl))
  
  for (i in 1:length(pred_mrl)){
    brms_b[i] = lt$t[which.min(abs(lt$mrl - pred_mrl[i]))]
  }
  
  # Calculate rmses
  pred_RMSE = sqrt(mean((df_sim_test$mrl - df_sim_test$pred_mrl)^2))
  beta_RMSE = sqrt(mean((betahat - true_betas$beta)^2))
  bioage_RMSE = sqrt(mean((brms_b - df_sim_test$b_age)^2))

  # MRL prediction plot
  plt = df_sim_test %>% ggplot(aes(x = mrl, y = pred_mrl)) +
    geom_point() +
    geom_abline(color = 'red') +
    labs(title = "True vs Predicted MRL for BRMS AFT Gompertz",
         subtitle = glue("RMSE = {pred_RMSE}")) +
    theme_minimal()
  # print(plt)

  # Beta comparison plot
  true_betas$betahat = betahat

  plt = true_betas %>%
    mutate(index = 1:p) %>%
    ggplot(aes(x = index)) +
    geom_segment(aes(xend = index, y = beta, yend = betahat), color = "gray") +
    geom_point(aes(y = beta, color = "True")) +
    geom_point(aes(y = betahat, color = "Predicted")) +
    scale_color_manual(values = c("True" = "black", "Predicted" = "skyblue"),
                       name = "Type") +
    ylab("β") +
    xlab("Index") +
    labs(title = "True vs Predicted β-coefficients for BRMS AFT Gompertz", 
         subtitle = glue("RMSE = {beta_RMSE %>% round(4)}")) +
    theme_minimal()
  # print(plt)
  
  return(list(
    test_data = df_sim_test,
    true_betas = true_betas,
    predicted_betas = betahat,
    predicted_bioage = brms_b,
    rmses = list(coef = beta_RMSE, mrl = pred_RMSE, bioage = bioage_RMSE),
    method = "BRMS gompertz",
    distr_ests = list(mu = mu_est, gamma = gamma_est, a = a_est, b = b_est)
  ))
}

3.3.2 Weibull implementation

brms_reg_weibull = function(df_sim, p, g, seed = 123) {
  true_betas = df_sim$true_betas
  lt = df_sim$lt
  df_sim = df_sim$df
  
  # Train/test split 50/50
  train_indices = sample(1:nrow(df_sim), nrow(df_sim)/2)
  df_sim_train = df_sim[train_indices,]
  df_sim_test = df_sim[-train_indices,]

  #### FIT MODEL
  preds = paste(paste0("x", 1:p), collapse = " + ")
  formula = as.formula(paste("age_end | cens(1 - status) + trunc(lb = age_start) ~", preds))

  priors = c(
    #prior(normal(0,0.5), class = "b"), # I took a wide sd due to convergence problems
    prior(horseshoe(df = 3, par_ratio = 0.1), class = "b"), # WOW!!!
    prior(normal(log(80), 0.2), class = Intercept), # about human life expectancy
    # prior(normal(10, 4), class = shape) # Hazard: bad, gamma better?
    prior(gamma(50, 5), class = shape) # not log scale!! (mean of gamma is a/b)
  )


  fit_bayes_aft = brm(
    formula = formula,
    data = df_sim_train,
    family = weibull(),
    prior = priors,
    # warmup = 250, 
    iter = 6000, chains = 2, cores = 2,
    control = list(adapt_delta = 0.97, max_treedepth = 15, stepsize = 0.01),
    seed = seed,
    # install.packages("cmdstanr", repos = c('https://stan-dev.r-universe.dev', getOption("repos")))
    backend = "cmdstanr",
    threads = threading(2),
    save_pars = save_pars(all = TRUE), init = 0

  )
  # Whenever you see the warning
  # "There were x divergent transitions after warmup.", you should really think about
  # increasing adapt_delta. To do this, write control = list(adapt_delta = <x>), where
  # <x> should usually be a value between 0.8 (current default) and 1. Increasing adapt_delta
  # will slow down the sampler but will decrease the number of divergent transitions threatening
  # the validity of your posterior samples.
  # https://cran.r-project.org/web/packages/brms/vignettes/brms_overview.pdf
  
  fb <<- fit_bayes_aft
  
  
  # Extract posterior samples and calculate estimates
  post = posterior_samples(fit_bayes_aft)
  beta_cols = paste0("b_x", 1:p)
  betahat = apply(post[, beta_cols], 2, mean)
  
  X_test = df_sim_test[,1:p] %>% as.matrix
  linpred = X_test %*% betahat
  
  shape_est = (mean(post$shape)) # Not exp because not in log scale????
  cat("Shape estimate:", shape_est, "\n")
  scale_est = exp(mean(post$b_Intercept))
  cat("Scale estimate", (scale_est), "\n")
  
  #### ESTIMATE MRL
  #  scale^(-shape) == lambda so we use it as est (from marijes script)

  lambda_est <- scale_est^(-shape_est)
  nu_est <- shape_est
  
  for (i in 1:nrow(df_sim_test)) {
    scale_adj = scale_est * exp(linpred[i])
    
    mrl_uncon <- integrate(function(t) {
      1 - pweibull(t, shape = shape_est, scale = scale_adj)
    }, lower = df_sim_test$age_start[i], upper = Inf)$value
    
    # Conditional survival at age_start
    s_cond <- 1 - pweibull(df_sim_test$age_start[i], shape = shape_est, scale = scale_adj)
    
    # Mean residual life
    df_sim_test$pred_mrl[i] <- mrl_uncon / s_cond

    # mrl_uncon <- integrate(weib_baseline_surv,
    #                        lower = (df_sim_test$age_start[i] * exp(linpred[i])), 
    #                      upper=Inf, lambda = lambda_est, nu = nu_est)$value
    # s_cond <-  weib_baseline_surv(df_sim_test$age_start[i] * exp(linpred[i]), 
    #                              lambda = lambda_est, nu = nu_est) 
    # df_sim_test$pred_mrl[i] <- (mrl_uncon / s_cond) * exp(-linpred[i])
  }

  # obtain biological age, assuming true lifetable is known
  pred_mrl = df_sim_test$pred_mrl
  brms_b = vector(length = length(pred_mrl))
  
  for (i in 1:length(pred_mrl)){
    brms_b[i] = lt$t[which.min(abs(lt$mrl - pred_mrl[i]))]
  }
  
  # Calculate rmses
  pred_RMSE = sqrt(mean((df_sim_test$mrl - df_sim_test$pred_mrl)^2))
  beta_RMSE = sqrt(mean((betahat - true_betas$beta)^2))
  bioage_RMSE = sqrt(mean((brms_b - df_sim_test$b_age)^2))

  # MRL prediction plot
  # plt = df_sim_test %>% ggplot(aes(x = mrl, y = pred_mrl)) +
  #   geom_point() +
  #   geom_abline(color = 'red') +
  #   labs(title = "True vs Predicted MRL for BRMS AFT Weibull",
  #        subtitle = glue("RMSE = {pred_RMSE}")) +
  #   theme_minimal()
  # print(plt)

  # Beta comparison plot
  true_betas$betahat = betahat

  # plt = true_betas %>%
  #   mutate(index = 1:p) %>%
  #   ggplot(aes(x = index)) +
  #   geom_segment(aes(xend = index, y = beta, yend = betahat), color = "gray") +
  #   geom_point(aes(y = beta, color = "True")) +
  #   geom_point(aes(y = betahat, color = "Predicted")) +
  #   scale_color_manual(values = c("True" = "black", "Predicted" = "skyblue"),
  #                      name = "Type") +
  #   ylab("β") +
  #   xlab("Index") +
  #   labs(title = "True vs Predicted β-coefficients for BRMS AFT Weibull",
  #        subtitle = glue("RMSE = {beta_RMSE %>% round(4)}")) +
  #   theme_minimal()
  # print(plt)
  
  return(list(
    test_data = df_sim_test,
    true_betas = true_betas,
    predicted_betas = betahat,
    predicted_bioage = brms_b,
    rmses = list(coef = beta_RMSE, mrl = pred_RMSE, bioage = bioage_RMSE),
    method = "BRMS weibull",
    distr_ests = list(shape = shape_est, scale = scale_est)
  ))
}

3.4 Monte Carlo Experiments

run_grid_experiments_mc = function(
  n_grid = c(200, 500, 1000),
  p_grid = c(20, 60, 100),
  g_grid = c(4, 6, 10),
  non_zero_grid = c(0.5, 0.75, 1.0),
  M = 10,
  pop_n = 100000,
  bsep = 4,
  X_rho = 0,
  target_snr = 0.1,
  method = "weibull",
  seed_base = 123
) {
  
  # Create parameter grid
  param_grid = expand.grid(
    n_obs = n_grid,
    p = p_grid,
    g = g_grid,
    non_zero = non_zero_grid
  )
  
  # Filter invalid combinations, p should be divisible by g (not irrational)
  param_grid = param_grid[param_grid$p %% param_grid$g == 0, ]
  
  cat(glue("Running {nrow(param_grid)} parameter combinations with {M} replications each\n"))
  cat(glue("Total experiments: {nrow(param_grid) * M}\n\n"))
  
  all_results = list()
  
  if (method == "gompertz") {
    a = exp(-9)
    b = 0.085
  } else if (method == "weibull") {
    a = 5 # Shape
    b = 80 # Scale
  }
  
  # Iterate through the grid
  for (i in 1:nrow(param_grid)) {
    # Params to compare
    params = param_grid[i,]
    
    cat(glue("\n<<<<<<<< Parameter Set {i}/{nrow(param_grid)} >>>>>>>>>>"), "\n")
    cat(glue("n = {params$n_obs} | p = {params$p} | g = {params$g} | non-zero = {params$non_zero}"), "\n")
    
    # Run M Monte Carlo replications
    mc_results = list()
    
    for (m in 1:M) {
      seed = seed_base + m
      
      cat(glue("  Rep {m}/{M}... "), "\n")
      
      # Generate data
      gres = generate_data(
        pop_n = pop_n, 
        obs_n = params$n_obs, 
        p = params$p,
        g = params$g,
        ltname = glue("mc_{.....}", "\n"),
        seed = seed, 
        force_recalc = FALSE, # IMPORTANT! SAME LT!
        X_plots = FALSE, 
        X_rho = X_rho, 
        β_rho = 0.9, 
        non_zero_betas = params$non_zero, 
        X_scale = 1,
        active_hazard = 0, 
        betasep = bsep, 
        bscale = 1, 
        betaplot = FALSE, 
        target_snr = target_snr,
        method = method,
        a = a,
        b = b
      )
      
      results = list()
      
      # AFT
      print("Running Frequentist AFT...")
      set.seed(seed)
      results$aft = tryCatch(
        aft_reg(df_sim = gres, p = params$p, method = method),
        error = function(e) NULL
      )
      
      # # PSBC
      # print("Running PSBC Bayesian AFT...")
      # set.seed(seed)
      # results$psbc = tryCatch(
      #   psbc_reg(df_sim = df_sim, p = params$p, g = params$g),
      #   error = function(e) NULL
      # )
      
      # BRMS
      print("Running BRMS Bayesian AFT...")
      set.seed(seed)
      results$brms = tryCatch(
        brms_reg(df_sim = gres, p = params$p, g = params$g, method = method),
        error = function(e) NULL
      )
      
      mc_results[[m]] = list(
        params = params,
        seed = seed,
        results = results
      )
      
      cat("Done\n")
    }
    ## TODO: SAVE RESULTS TO FILES IN SOME EXPERIMENT FOLDER
    all_results[[i]] = list(
      params = params,
      mc_results = mc_results
    )
  }
  
  return(all_results)
}

4 Experiments Monte Carlo

options(mc.cores = 4)

4.1 \(n \in \{200, 400, 1000\}, p \in \{20\}, g = 5\)

# results = run_grid_experiments_mc(n_grid = c(200), p_grid = c(10), g_grid = c(5), non_zero_grid = c(1), M = 1)
# rsp = parse_results(all_results = results)
# plot_results_by_n(all_results = results, target_p = 20)

4.2 \(n \in \{200, 400, 1000\}, p \in \{200\}, g = 5\)

# results = run_grid_experiments_mc(n_grid = c(400, 800, 2000), p_grid = c(200), g_grid = c(5), non_zero_grid = c(1), M = 1)
# plot_results_by_n(all_results = results, target_p = 200)

5 High dimensional example (p = n)

Since obs_n = 200, train test split of 50/50 means that it is trained on 100 observations.

# rmses = run_experiment(p = 10, g = 4, ltname = "p60_may20", seed = 43, nzb = 0.75, bsep = 4, X_rho = 0, obs_n = 200)
# print(rmses)

Small tests:

6 DEBUG - before performing actual MC experiments

6.1 Debug function:

debug_methods = function(n = 100, p = 10, g = 5, method = "weibull", seed = 123, t_SNR = 2, perf_A = T, nzg = 1) {
  
  cat("**** SETTINGS ****\n")
  cat("n =", n, "p =", p, "g =", g, "method =", method, "\n\n")
  
  cat("Generating data...\n")
  if (method == "weibull") {
    a = 10 # shape
    b = 80 # scale
  } else {
    a = exp(-9) # gompertz a
    b = 0.085 # gompertz b
  }
  
  gres = generate_data(
    pop_n = as.integer(1e4), # 10,000
    obs_n = n * 2, # x2 due to train test split
    p = p,
    g = g,
    ltname = glue("debug_{method}"),
    seed = seed,
    force_recalc = TRUE,
    X_plots = FALSE,
    X_rho = 0,
    β_rho = 0.9,
    non_zero_betas = nzg,
    X_scale = 1,
    active_hazard = 0,
    betasep = 4,
    bscale = 1,
    betaplot = TRUE,
    target_snr = t_SNR,
    method = method,
    a = a,
    b = b
  )
  
  cat("Data generated successfully!\n")
  cat("Dataset size:", nrow(gres$df), "x", ncol(gres$df), "\n\n")
  
  if (!perf_A) {
    return(gres)
  }
  
  cat("**** TESTING AFT ****\n")
  set.seed(seed)
  aft_result = aft_reg(df_sim = gres, p = p, method = method)
  
  if (!is.null(aft_result)) {
    cat("AFT SUCCESS!\n")
    cat("AFT RMSEs:", "\n")
    print(aft_result$rmses)
    plot(fa)
  } else {
    cat("AFT FAILED!\n")
  }
  # aft_result = NULL # DEBUG
  
  cat("\n**** TESTING BRMS ****\n")
  set.seed(seed)
  brms_result = brms_reg(df_sim = gres, p = p, g = g, method = method, seed = seed)
  
  if (!is.null(brms_result)) {
    cat("BRMS SUCCESS!\n") 
    cat("BRMS RMSEs:", "\n")
    print(brms_result$rmses)
    plot(fb, variable = c("b_Intercept", "shape", "b_x1"))
  } else {
    cat("BRMS FAILED!\n")
  }
  
  res = list(aft = aft_result, brms = brms_result)
  cat("\n**** COMPARISON ****\n")
  if (!is.null(brms_result)) {
    if (is.null(aft_result)) {
      res$aft = list(
        test_data = data.frame(
          mrl = rep(0, n),
          pred_mrl = rep(0, n),
          b_age = rep(0, n)
        ),
        true_betas = res$brms$true_betas$beta,
        predicted_betas = rep(0, p),
        predicted_bioage = rep(0, n),
        rmses = list(coef = 0, mrl = 0, bioage = 0),
        method = "AFT (failed)"
      )
      aft_result = res$aft
    }
    
    p_vars = length(res$brms$true_betas$beta)
    beta_comparison = data.frame(
      index = 1:p_vars,
      true_beta = res$brms$true_betas$beta,
      aft_beta = res$aft$predicted_betas,
      brms_beta = res$brms$predicted_betas
    )
    
    # Beta plot
    p1 = ggplot(beta_comparison, aes(x = index)) +
      geom_segment(aes(xend = index, y = true_beta, yend = aft_beta), 
                   color = "lightgray", alpha = 0.5) +
      geom_segment(aes(xend = index, y = true_beta, yend = brms_beta), 
                   color = "lightblue", alpha = 0.5) +
      geom_point(aes(y = true_beta, color = "True"), size = 1.5) +
      geom_point(aes(y = aft_beta, color = "AFT"), size = 1) +
      geom_point(aes(y = brms_beta, color = "BRMS"), size = 1) +
      scale_color_manual(values = c("True" = "black", "AFT" = "red", "BRMS" = "blue")) +
      geom_hline(yintercept = 0) +
      labs(title = "Beta Coefficient Comparison: AFT vs BRMS",
           subtitle = paste("AFT Beta RMSE:", round(res$aft$rmses$coef, 4), 
                           "| BRMS Beta RMSE:", round(res$brms$rmses$coef, 4)),
           x = "Coefficient Index", y = "Beta Value") +
      theme_classic()
    
    print(p1)
    
    n_test = min(nrow(res$aft$test_data), nrow(res$brms$test_data))
    
    mrl_comparison = data.frame(
      true_mrl = res$brms$test_data$mrl[1:n_test],
      aft_mrl = res$aft$test_data$pred_mrl[1:n_test],
      brms_mrl = res$brms$test_data$pred_mrl[1:n_test]
    )
    
    mrl_comparison = mrl_comparison[complete.cases(mrl_comparison), ]
    
    p2 = ggplot(mrl_comparison, aes(x = true_mrl)) +
      geom_point(aes(y = aft_mrl, color = "AFT"), alpha = 0.6) +
      geom_point(aes(y = brms_mrl, color = "BRMS"), alpha = 0.6) +
      geom_abline(intercept = 0, slope = 1, color = "black", linetype = "dashed") +
      scale_color_manual(values = c("AFT" = "red", "BRMS" = "blue")) +
      labs(title = "MRL Prediction Comparison: AFT vs BRMS",
           x = "True MRL", y = "Predicted MRL",
           subtitle = paste("AFT RMSE:", round(res$aft$rmses$mrl, 3), 
                           "| BRMS RMSE:", round(res$brms$rmses$mrl, 3))) +
      theme_minimal()
    
    print(p2)
    
    # Bioage
    bioage_comparison = data.frame(
      true_b_age = brms_result$test_data$b_age,
      aft_b_age = aft_result$predicted_bioage,
      brms_b_age = brms_result$predicted_bioage
    )
    
    bioage_comparison = bioage_comparison[complete.cases(bioage_comparison), ]
    
    p3 = ggplot(bioage_comparison, aes(x = true_b_age)) +
      geom_point(aes(y = aft_b_age, color = "AFT"), alpha = 0.6) +
      geom_point(aes(y = brms_b_age, color = "BRMS"), alpha = 0.6) +
      geom_abline(intercept = 0, slope = 1, color = "black", linetype = "dashed") +
      scale_color_manual(values = c("AFT" = "red", "BRMS" = "blue")) +
      labs(title = "Biological Age Prediction Comparison: AFT vs BRMS",
           x = "True Biological Age", y = "Predicted Biological Age",
           subtitle = paste("AFT RMSE:", round(aft_result$rmses$bioage, 3), 
                           "| BRMS RMSE:", round(brms_result$rmses$bioage, 3))) +
      theme_minimal()
    
    print(p3)
        
    cat("\n**** SUMMARY COMPARISON ****\n")
    cat("Beta RMSE - AFT:", round(res$aft$rmses$coef, 4), "\n")
    cat("Beta RMSE - BRMS:", round(res$brms$rmses$coef, 4), "\n")
    cat("MRL RMSE - AFT:", round(res$aft$rmses$mrl, 4), "\n") 
    cat("MRL RMSE - BRMS:", round(res$brms$rmses$mrl, 4), "\n")
    cat("BioAge RMSE - AFT", round(res$aft$rmses$bioage, 4), "\n")
    cat("BioAge RMSE - BRMS", round(res$brms$rmses$bioage, 4), "\n")
  }
  
  return(res)
}

6.2 Method: Weibull; p = 20; \(n \in \{250, 500\}\); \(r_{nz} \in \{ 1\}\)

6.2.1 n = 250

source("R/generate_data.R")

debug_methods(n = 250, p = 20, g = 5, method = "weibull", seed = 40, t_SNR = 4, perf_A = T, nzg = 1) %>%
  system.time()
## **** SETTINGS ****
## n = 250 p = 20 g = 5 method = weibull 
## 
## Generating data...
## [1] "Generating βs..."
## Mean of betas pre-scaling: 0.286206588904146
## Lower range of beta values pre-scaling: -3.42901337506225
## Upper range of beta values pre-scaling: 4.47580730241499
## Warning: `as.tibble()` was deprecated in tibble 2.0.0.
## ℹ Please use `as_tibble()` instead.
## ℹ The signature and semantics have changed, see `?as_tibble`.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if
## `.name_repair` is omitted as of tibble 2.0.0.
## ℹ Using compatibility `.name_repair`.
## ℹ The deprecated feature was likely used in the tibble package.
##   Please report the issue at <https://github.com/tidyverse/tibble/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Method: weibull
## Target SNR: 4
## Epsilon variance: 0.0165
## Mean of betas post-scaling: 2.25514051876985e-18
## Lower range of beta values post-scaling: -0.0422802473570912
## Upper range of beta values post-scaling: 0.0476788336122017
## Achieved SNR: 4

## [1] "Generating X..."
## [1] "Done!"
## [1] "Generating Weibull population lifetable..."
## Range of death: 27.7661009696673
## Range of death: 143.007815466488

## [1] "Done!"

## [1] "Generating data from Weibull lifetable..."
## Nr of valid indices: 2179
## Range of linpred values: -0.5733417 0.455935 
## [1] "Done!"
## Data generated successfully!
## Dataset size: 500 x 29 
## 
## **** TESTING AFT ****
## [1] "----------"
## shape_est: 11.90365 
## scale_est: 78.61759 
## [1] "----------"
## AFT SUCCESS!
## AFT RMSEs: 
## $coef
## [1] 0.01134541
## 
## $mrl
## [1] 3.751361
## 
## $bioage
## [1] 6.265548

## 
## **** TESTING BRMS ****
## Start sampling
## Warning in .fun(data = .x1, seed = .x2, init = .x3, adapt_delta = .x4,
## max_treedepth = .x5, : 'stepsize' is deprecated. Please use 'step_size'
## instead.
## Running MCMC with 2 parallel chains, with 2 thread(s) per chain...
## 
## Chain 1 Iteration:    1 / 6000 [  0%]  (Warmup) 
## Chain 2 Iteration:    1 / 6000 [  0%]  (Warmup) 
## Chain 2 Iteration:  100 / 6000 [  1%]  (Warmup) 
## Chain 1 Iteration:  100 / 6000 [  1%]  (Warmup) 
## Chain 2 Iteration:  200 / 6000 [  3%]  (Warmup) 
## Chain 2 Iteration:  300 / 6000 [  5%]  (Warmup) 
## Chain 1 Iteration:  200 / 6000 [  3%]  (Warmup) 
## Chain 1 Iteration:  300 / 6000 [  5%]  (Warmup) 
## Chain 2 Iteration:  400 / 6000 [  6%]  (Warmup) 
## Chain 1 Iteration:  400 / 6000 [  6%]  (Warmup) 
## Chain 2 Iteration:  500 / 6000 [  8%]  (Warmup) 
## Chain 1 Iteration:  500 / 6000 [  8%]  (Warmup) 
## Chain 2 Iteration:  600 / 6000 [ 10%]  (Warmup) 
## Chain 1 Iteration:  600 / 6000 [ 10%]  (Warmup) 
## Chain 2 Iteration:  700 / 6000 [ 11%]  (Warmup) 
## Chain 1 Iteration:  700 / 6000 [ 11%]  (Warmup) 
## Chain 1 Iteration:  800 / 6000 [ 13%]  (Warmup) 
## Chain 2 Iteration:  800 / 6000 [ 13%]  (Warmup) 
## Chain 1 Iteration:  900 / 6000 [ 15%]  (Warmup) 
## Chain 2 Iteration:  900 / 6000 [ 15%]  (Warmup) 
## Chain 1 Iteration: 1000 / 6000 [ 16%]  (Warmup) 
## Chain 2 Iteration: 1000 / 6000 [ 16%]  (Warmup) 
## Chain 1 Iteration: 1100 / 6000 [ 18%]  (Warmup) 
## Chain 2 Iteration: 1100 / 6000 [ 18%]  (Warmup) 
## Chain 1 Iteration: 1200 / 6000 [ 20%]  (Warmup) 
## Chain 2 Iteration: 1200 / 6000 [ 20%]  (Warmup) 
## Chain 1 Iteration: 1300 / 6000 [ 21%]  (Warmup) 
## Chain 2 Iteration: 1300 / 6000 [ 21%]  (Warmup) 
## Chain 1 Iteration: 1400 / 6000 [ 23%]  (Warmup) 
## Chain 2 Iteration: 1400 / 6000 [ 23%]  (Warmup) 
## Chain 1 Iteration: 1500 / 6000 [ 25%]  (Warmup) 
## Chain 2 Iteration: 1500 / 6000 [ 25%]  (Warmup) 
## Chain 1 Iteration: 1600 / 6000 [ 26%]  (Warmup) 
## Chain 2 Iteration: 1600 / 6000 [ 26%]  (Warmup) 
## Chain 1 Iteration: 1700 / 6000 [ 28%]  (Warmup) 
## Chain 2 Iteration: 1700 / 6000 [ 28%]  (Warmup) 
## Chain 1 Iteration: 1800 / 6000 [ 30%]  (Warmup) 
## Chain 2 Iteration: 1800 / 6000 [ 30%]  (Warmup) 
## Chain 1 Iteration: 1900 / 6000 [ 31%]  (Warmup) 
## Chain 2 Iteration: 1900 / 6000 [ 31%]  (Warmup) 
## Chain 1 Iteration: 2000 / 6000 [ 33%]  (Warmup) 
## Chain 2 Iteration: 2000 / 6000 [ 33%]  (Warmup) 
## Chain 1 Iteration: 2100 / 6000 [ 35%]  (Warmup) 
## Chain 2 Iteration: 2100 / 6000 [ 35%]  (Warmup) 
## Chain 1 Iteration: 2200 / 6000 [ 36%]  (Warmup) 
## Chain 2 Iteration: 2200 / 6000 [ 36%]  (Warmup) 
## Chain 1 Iteration: 2300 / 6000 [ 38%]  (Warmup) 
## Chain 2 Iteration: 2300 / 6000 [ 38%]  (Warmup) 
## Chain 1 Iteration: 2400 / 6000 [ 40%]  (Warmup) 
## Chain 2 Iteration: 2400 / 6000 [ 40%]  (Warmup) 
## Chain 1 Iteration: 2500 / 6000 [ 41%]  (Warmup) 
## Chain 2 Iteration: 2500 / 6000 [ 41%]  (Warmup) 
## Chain 1 Iteration: 2600 / 6000 [ 43%]  (Warmup) 
## Chain 2 Iteration: 2600 / 6000 [ 43%]  (Warmup) 
## Chain 1 Iteration: 2700 / 6000 [ 45%]  (Warmup) 
## Chain 2 Iteration: 2700 / 6000 [ 45%]  (Warmup) 
## Chain 1 Iteration: 2800 / 6000 [ 46%]  (Warmup) 
## Chain 2 Iteration: 2800 / 6000 [ 46%]  (Warmup) 
## Chain 2 Iteration: 2900 / 6000 [ 48%]  (Warmup) 
## Chain 1 Iteration: 2900 / 6000 [ 48%]  (Warmup) 
## Chain 1 Iteration: 3000 / 6000 [ 50%]  (Warmup) 
## Chain 1 Iteration: 3001 / 6000 [ 50%]  (Sampling) 
## Chain 2 Iteration: 3000 / 6000 [ 50%]  (Warmup) 
## Chain 2 Iteration: 3001 / 6000 [ 50%]  (Sampling) 
## Chain 1 Iteration: 3100 / 6000 [ 51%]  (Sampling) 
## Chain 2 Iteration: 3100 / 6000 [ 51%]  (Sampling) 
## Chain 1 Iteration: 3200 / 6000 [ 53%]  (Sampling) 
## Chain 2 Iteration: 3200 / 6000 [ 53%]  (Sampling) 
## Chain 1 Iteration: 3300 / 6000 [ 55%]  (Sampling) 
## Chain 2 Iteration: 3300 / 6000 [ 55%]  (Sampling) 
## Chain 1 Iteration: 3400 / 6000 [ 56%]  (Sampling) 
## Chain 2 Iteration: 3400 / 6000 [ 56%]  (Sampling) 
## Chain 1 Iteration: 3500 / 6000 [ 58%]  (Sampling) 
## Chain 2 Iteration: 3500 / 6000 [ 58%]  (Sampling) 
## Chain 1 Iteration: 3600 / 6000 [ 60%]  (Sampling) 
## Chain 2 Iteration: 3600 / 6000 [ 60%]  (Sampling) 
## Chain 1 Iteration: 3700 / 6000 [ 61%]  (Sampling) 
## Chain 2 Iteration: 3700 / 6000 [ 61%]  (Sampling) 
## Chain 1 Iteration: 3800 / 6000 [ 63%]  (Sampling) 
## Chain 2 Iteration: 3800 / 6000 [ 63%]  (Sampling) 
## Chain 1 Iteration: 3900 / 6000 [ 65%]  (Sampling) 
## Chain 2 Iteration: 3900 / 6000 [ 65%]  (Sampling) 
## Chain 1 Iteration: 4000 / 6000 [ 66%]  (Sampling) 
## Chain 2 Iteration: 4000 / 6000 [ 66%]  (Sampling) 
## Chain 1 Iteration: 4100 / 6000 [ 68%]  (Sampling) 
## Chain 2 Iteration: 4100 / 6000 [ 68%]  (Sampling) 
## Chain 1 Iteration: 4200 / 6000 [ 70%]  (Sampling) 
## Chain 2 Iteration: 4200 / 6000 [ 70%]  (Sampling) 
## Chain 1 Iteration: 4300 / 6000 [ 71%]  (Sampling) 
## Chain 2 Iteration: 4300 / 6000 [ 71%]  (Sampling) 
## Chain 1 Iteration: 4400 / 6000 [ 73%]  (Sampling) 
## Chain 2 Iteration: 4400 / 6000 [ 73%]  (Sampling) 
## Chain 1 Iteration: 4500 / 6000 [ 75%]  (Sampling) 
## Chain 2 Iteration: 4500 / 6000 [ 75%]  (Sampling) 
## Chain 1 Iteration: 4600 / 6000 [ 76%]  (Sampling) 
## Chain 2 Iteration: 4600 / 6000 [ 76%]  (Sampling) 
## Chain 1 Iteration: 4700 / 6000 [ 78%]  (Sampling) 
## Chain 2 Iteration: 4700 / 6000 [ 78%]  (Sampling) 
## Chain 1 Iteration: 4800 / 6000 [ 80%]  (Sampling) 
## Chain 2 Iteration: 4800 / 6000 [ 80%]  (Sampling) 
## Chain 1 Iteration: 4900 / 6000 [ 81%]  (Sampling) 
## Chain 2 Iteration: 4900 / 6000 [ 81%]  (Sampling) 
## Chain 1 Iteration: 5000 / 6000 [ 83%]  (Sampling) 
## Chain 2 Iteration: 5000 / 6000 [ 83%]  (Sampling) 
## Chain 1 Iteration: 5100 / 6000 [ 85%]  (Sampling) 
## Chain 2 Iteration: 5100 / 6000 [ 85%]  (Sampling) 
## Chain 1 Iteration: 5200 / 6000 [ 86%]  (Sampling) 
## Chain 2 Iteration: 5200 / 6000 [ 86%]  (Sampling) 
## Chain 1 Iteration: 5300 / 6000 [ 88%]  (Sampling) 
## Chain 2 Iteration: 5300 / 6000 [ 88%]  (Sampling) 
## Chain 1 Iteration: 5400 / 6000 [ 90%]  (Sampling) 
## Chain 2 Iteration: 5400 / 6000 [ 90%]  (Sampling) 
## Chain 1 Iteration: 5500 / 6000 [ 91%]  (Sampling) 
## Chain 2 Iteration: 5500 / 6000 [ 91%]  (Sampling) 
## Chain 1 Iteration: 5600 / 6000 [ 93%]  (Sampling) 
## Chain 2 Iteration: 5600 / 6000 [ 93%]  (Sampling) 
## Chain 1 Iteration: 5700 / 6000 [ 95%]  (Sampling) 
## Chain 2 Iteration: 5700 / 6000 [ 95%]  (Sampling) 
## Chain 1 Iteration: 5800 / 6000 [ 96%]  (Sampling) 
## Chain 2 Iteration: 5800 / 6000 [ 96%]  (Sampling) 
## Chain 1 Iteration: 5900 / 6000 [ 98%]  (Sampling) 
## Chain 2 Iteration: 5900 / 6000 [ 98%]  (Sampling) 
## Chain 1 Iteration: 6000 / 6000 [100%]  (Sampling) 
## Chain 1 finished in 31.8 seconds.
## Chain 2 Iteration: 6000 / 6000 [100%]  (Sampling) 
## Chain 2 finished in 31.9 seconds.
## 
## Both chains finished successfully.
## Mean chain execution time: 31.9 seconds.
## Total execution time: 32.0 seconds.
## Warning: 19 of 6000 (0.0%) transitions ended with a divergence.
## See https://mc-stan.org/misc/warnings for details.
## Loading required package: rstan
## Warning: package 'rstan' was built under R version 4.3.3
## Loading required package: StanHeaders
## Warning: package 'StanHeaders' was built under R version 4.3.3
## 
## rstan version 2.32.7 (Stan version 2.32.2)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
## For within-chain threading using `reduce_sum()` or `map_rect()` Stan functions,
## change `threads_per_chain` option:
## rstan_options(threads_per_chain = 1)
## 
## Attaching package: 'rstan'
## The following object is masked from 'package:tidyr':
## 
##     extract
## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.
## Shape estimate: 9.9499 
## Scale estimate 75.63284 
## BRMS SUCCESS!
## BRMS RMSEs: 
## $coef
## [1] 0.01258627
## 
## $mrl
## [1] 5.790734
## 
## $bioage
## [1] 8.876742

## 
## **** COMPARISON ****

## 
## **** SUMMARY COMPARISON ****
## Beta RMSE - AFT: 0.0113 
## Beta RMSE - BRMS: 0.0126 
## MRL RMSE - AFT: 3.7514 
## MRL RMSE - BRMS: 5.7907 
## BioAge RMSE - AFT 6.2655 
## BioAge RMSE - BRMS 8.8767
##    user  system elapsed 
## 117.908  11.205  45.829

6.2.2 n = 500

debug_methods(n = 500, p = 20, g = 5, method = "weibull", seed = 40, t_SNR = 4, perf_A = T, nzg = 1) %>%
  system.time()
## **** SETTINGS ****
## n = 500 p = 20 g = 5 method = weibull 
## 
## Generating data...
## [1] "Generating βs..."
## Mean of betas pre-scaling: 0.286206588904146
## Lower range of beta values pre-scaling: -3.42901337506225
## Upper range of beta values pre-scaling: 4.47580730241499
## Method: weibull
## Target SNR: 4
## Epsilon variance: 0.0165
## Mean of betas post-scaling: 2.25514051876985e-18
## Lower range of beta values post-scaling: -0.0422802473570912
## Upper range of beta values post-scaling: 0.0476788336122017
## Achieved SNR: 4

## [1] "Generating X..."
## [1] "Done!"
## [1] "Generating Weibull population lifetable..."
## Range of death: 27.7661009696673
## Range of death: 143.007815466488

## [1] "Done!"

## [1] "Generating data from Weibull lifetable..."
## Nr of valid indices: 4407
## Range of linpred values: -0.5280914 0.5125125 
## [1] "Done!"
## Data generated successfully!
## Dataset size: 1000 x 29 
## 
## **** TESTING AFT ****
## [1] "----------"
## shape_est: 10.06515 
## scale_est: 80.82797 
## [1] "----------"
## AFT SUCCESS!
## AFT RMSEs: 
## $coef
## [1] 0.009912031
## 
## $mrl
## [1] 3.294007
## 
## $bioage
## [1] 6.440819

## 
## **** TESTING BRMS ****
## Start sampling
## Warning in .fun(data = .x1, seed = .x2, init = .x3, adapt_delta = .x4,
## max_treedepth = .x5, : 'stepsize' is deprecated. Please use 'step_size'
## instead.
## Running MCMC with 2 parallel chains, with 2 thread(s) per chain...
## 
## Chain 1 Iteration:    1 / 6000 [  0%]  (Warmup) 
## Chain 2 Iteration:    1 / 6000 [  0%]  (Warmup) 
## Chain 1 Iteration:  100 / 6000 [  1%]  (Warmup) 
## Chain 2 Iteration:  100 / 6000 [  1%]  (Warmup) 
## Chain 1 Iteration:  200 / 6000 [  3%]  (Warmup) 
## Chain 2 Iteration:  200 / 6000 [  3%]  (Warmup) 
## Chain 1 Iteration:  300 / 6000 [  5%]  (Warmup) 
## Chain 2 Iteration:  300 / 6000 [  5%]  (Warmup) 
## Chain 1 Iteration:  400 / 6000 [  6%]  (Warmup) 
## Chain 2 Iteration:  400 / 6000 [  6%]  (Warmup) 
## Chain 2 Iteration:  500 / 6000 [  8%]  (Warmup) 
## Chain 1 Iteration:  500 / 6000 [  8%]  (Warmup) 
## Chain 1 Iteration:  600 / 6000 [ 10%]  (Warmup) 
## Chain 2 Iteration:  600 / 6000 [ 10%]  (Warmup) 
## Chain 2 Iteration:  700 / 6000 [ 11%]  (Warmup) 
## Chain 1 Iteration:  700 / 6000 [ 11%]  (Warmup) 
## Chain 2 Iteration:  800 / 6000 [ 13%]  (Warmup) 
## Chain 1 Iteration:  800 / 6000 [ 13%]  (Warmup) 
## Chain 2 Iteration:  900 / 6000 [ 15%]  (Warmup) 
## Chain 1 Iteration:  900 / 6000 [ 15%]  (Warmup) 
## Chain 2 Iteration: 1000 / 6000 [ 16%]  (Warmup) 
## Chain 1 Iteration: 1000 / 6000 [ 16%]  (Warmup) 
## Chain 2 Iteration: 1100 / 6000 [ 18%]  (Warmup) 
## Chain 1 Iteration: 1100 / 6000 [ 18%]  (Warmup) 
## Chain 2 Iteration: 1200 / 6000 [ 20%]  (Warmup) 
## Chain 1 Iteration: 1200 / 6000 [ 20%]  (Warmup) 
## Chain 2 Iteration: 1300 / 6000 [ 21%]  (Warmup) 
## Chain 1 Iteration: 1300 / 6000 [ 21%]  (Warmup) 
## Chain 2 Iteration: 1400 / 6000 [ 23%]  (Warmup) 
## Chain 1 Iteration: 1400 / 6000 [ 23%]  (Warmup) 
## Chain 2 Iteration: 1500 / 6000 [ 25%]  (Warmup) 
## Chain 1 Iteration: 1500 / 6000 [ 25%]  (Warmup) 
## Chain 2 Iteration: 1600 / 6000 [ 26%]  (Warmup) 
## Chain 1 Iteration: 1600 / 6000 [ 26%]  (Warmup) 
## Chain 2 Iteration: 1700 / 6000 [ 28%]  (Warmup) 
## Chain 1 Iteration: 1700 / 6000 [ 28%]  (Warmup) 
## Chain 2 Iteration: 1800 / 6000 [ 30%]  (Warmup) 
## Chain 1 Iteration: 1800 / 6000 [ 30%]  (Warmup) 
## Chain 2 Iteration: 1900 / 6000 [ 31%]  (Warmup) 
## Chain 1 Iteration: 1900 / 6000 [ 31%]  (Warmup) 
## Chain 2 Iteration: 2000 / 6000 [ 33%]  (Warmup) 
## Chain 1 Iteration: 2000 / 6000 [ 33%]  (Warmup) 
## Chain 2 Iteration: 2100 / 6000 [ 35%]  (Warmup) 
## Chain 1 Iteration: 2100 / 6000 [ 35%]  (Warmup) 
## Chain 2 Iteration: 2200 / 6000 [ 36%]  (Warmup) 
## Chain 1 Iteration: 2200 / 6000 [ 36%]  (Warmup) 
## Chain 2 Iteration: 2300 / 6000 [ 38%]  (Warmup) 
## Chain 1 Iteration: 2300 / 6000 [ 38%]  (Warmup) 
## Chain 2 Iteration: 2400 / 6000 [ 40%]  (Warmup) 
## Chain 1 Iteration: 2400 / 6000 [ 40%]  (Warmup) 
## Chain 2 Iteration: 2500 / 6000 [ 41%]  (Warmup) 
## Chain 1 Iteration: 2500 / 6000 [ 41%]  (Warmup) 
## Chain 2 Iteration: 2600 / 6000 [ 43%]  (Warmup) 
## Chain 1 Iteration: 2600 / 6000 [ 43%]  (Warmup) 
## Chain 2 Iteration: 2700 / 6000 [ 45%]  (Warmup) 
## Chain 1 Iteration: 2700 / 6000 [ 45%]  (Warmup) 
## Chain 2 Iteration: 2800 / 6000 [ 46%]  (Warmup) 
## Chain 1 Iteration: 2800 / 6000 [ 46%]  (Warmup) 
## Chain 2 Iteration: 2900 / 6000 [ 48%]  (Warmup) 
## Chain 1 Iteration: 2900 / 6000 [ 48%]  (Warmup) 
## Chain 2 Iteration: 3000 / 6000 [ 50%]  (Warmup) 
## Chain 2 Iteration: 3001 / 6000 [ 50%]  (Sampling) 
## Chain 1 Iteration: 3000 / 6000 [ 50%]  (Warmup) 
## Chain 1 Iteration: 3001 / 6000 [ 50%]  (Sampling) 
## Chain 2 Iteration: 3100 / 6000 [ 51%]  (Sampling) 
## Chain 1 Iteration: 3100 / 6000 [ 51%]  (Sampling) 
## Chain 2 Iteration: 3200 / 6000 [ 53%]  (Sampling) 
## Chain 1 Iteration: 3200 / 6000 [ 53%]  (Sampling) 
## Chain 2 Iteration: 3300 / 6000 [ 55%]  (Sampling) 
## Chain 1 Iteration: 3300 / 6000 [ 55%]  (Sampling) 
## Chain 2 Iteration: 3400 / 6000 [ 56%]  (Sampling) 
## Chain 1 Iteration: 3400 / 6000 [ 56%]  (Sampling) 
## Chain 2 Iteration: 3500 / 6000 [ 58%]  (Sampling) 
## Chain 1 Iteration: 3500 / 6000 [ 58%]  (Sampling) 
## Chain 2 Iteration: 3600 / 6000 [ 60%]  (Sampling) 
## Chain 1 Iteration: 3600 / 6000 [ 60%]  (Sampling) 
## Chain 2 Iteration: 3700 / 6000 [ 61%]  (Sampling) 
## Chain 1 Iteration: 3700 / 6000 [ 61%]  (Sampling) 
## Chain 2 Iteration: 3800 / 6000 [ 63%]  (Sampling) 
## Chain 1 Iteration: 3800 / 6000 [ 63%]  (Sampling) 
## Chain 2 Iteration: 3900 / 6000 [ 65%]  (Sampling) 
## Chain 1 Iteration: 3900 / 6000 [ 65%]  (Sampling) 
## Chain 2 Iteration: 4000 / 6000 [ 66%]  (Sampling) 
## Chain 1 Iteration: 4000 / 6000 [ 66%]  (Sampling) 
## Chain 2 Iteration: 4100 / 6000 [ 68%]  (Sampling) 
## Chain 2 Iteration: 4200 / 6000 [ 70%]  (Sampling) 
## Chain 1 Iteration: 4100 / 6000 [ 68%]  (Sampling) 
## Chain 2 Iteration: 4300 / 6000 [ 71%]  (Sampling) 
## Chain 1 Iteration: 4200 / 6000 [ 70%]  (Sampling) 
## Chain 2 Iteration: 4400 / 6000 [ 73%]  (Sampling) 
## Chain 1 Iteration: 4300 / 6000 [ 71%]  (Sampling) 
## Chain 2 Iteration: 4500 / 6000 [ 75%]  (Sampling) 
## Chain 1 Iteration: 4400 / 6000 [ 73%]  (Sampling) 
## Chain 2 Iteration: 4600 / 6000 [ 76%]  (Sampling) 
## Chain 1 Iteration: 4500 / 6000 [ 75%]  (Sampling) 
## Chain 2 Iteration: 4700 / 6000 [ 78%]  (Sampling) 
## Chain 1 Iteration: 4600 / 6000 [ 76%]  (Sampling) 
## Chain 2 Iteration: 4800 / 6000 [ 80%]  (Sampling) 
## Chain 1 Iteration: 4700 / 6000 [ 78%]  (Sampling) 
## Chain 2 Iteration: 4900 / 6000 [ 81%]  (Sampling) 
## Chain 1 Iteration: 4800 / 6000 [ 80%]  (Sampling) 
## Chain 2 Iteration: 5000 / 6000 [ 83%]  (Sampling) 
## Chain 1 Iteration: 4900 / 6000 [ 81%]  (Sampling) 
## Chain 2 Iteration: 5100 / 6000 [ 85%]  (Sampling) 
## Chain 1 Iteration: 5000 / 6000 [ 83%]  (Sampling) 
## Chain 2 Iteration: 5200 / 6000 [ 86%]  (Sampling) 
## Chain 1 Iteration: 5100 / 6000 [ 85%]  (Sampling) 
## Chain 2 Iteration: 5300 / 6000 [ 88%]  (Sampling) 
## Chain 1 Iteration: 5200 / 6000 [ 86%]  (Sampling) 
## Chain 2 Iteration: 5400 / 6000 [ 90%]  (Sampling) 
## Chain 1 Iteration: 5300 / 6000 [ 88%]  (Sampling) 
## Chain 2 Iteration: 5500 / 6000 [ 91%]  (Sampling) 
## Chain 1 Iteration: 5400 / 6000 [ 90%]  (Sampling) 
## Chain 2 Iteration: 5600 / 6000 [ 93%]  (Sampling) 
## Chain 1 Iteration: 5500 / 6000 [ 91%]  (Sampling) 
## Chain 2 Iteration: 5700 / 6000 [ 95%]  (Sampling) 
## Chain 1 Iteration: 5600 / 6000 [ 93%]  (Sampling) 
## Chain 2 Iteration: 5800 / 6000 [ 96%]  (Sampling) 
## Chain 1 Iteration: 5700 / 6000 [ 95%]  (Sampling) 
## Chain 2 Iteration: 5900 / 6000 [ 98%]  (Sampling) 
## Chain 1 Iteration: 5800 / 6000 [ 96%]  (Sampling) 
## Chain 2 Iteration: 6000 / 6000 [100%]  (Sampling) 
## Chain 2 finished in 79.4 seconds.
## Chain 1 Iteration: 5900 / 6000 [ 98%]  (Sampling) 
## Chain 1 Iteration: 6000 / 6000 [100%]  (Sampling) 
## Chain 1 finished in 80.8 seconds.
## 
## Both chains finished successfully.
## Mean chain execution time: 80.1 seconds.
## Total execution time: 80.9 seconds.
## Warning: 42 of 6000 (1.0%) transitions ended with a divergence.
## See https://mc-stan.org/misc/warnings for details.
## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.
## Shape estimate: 9.378186 
## Scale estimate 76.9985 
## BRMS SUCCESS!
## BRMS RMSEs: 
## $coef
## [1] 0.01035903
## 
## $mrl
## [1] 4.586242
## 
## $bioage
## [1] 8.699857

## 
## **** COMPARISON ****

## 
## **** SUMMARY COMPARISON ****
## Beta RMSE - AFT: 0.0099 
## Beta RMSE - BRMS: 0.0104 
## MRL RMSE - AFT: 3.294 
## MRL RMSE - BRMS: 4.5862 
## BioAge RMSE - AFT 6.4408 
## BioAge RMSE - BRMS 8.6999
##    user  system elapsed 
## 264.076  22.496  85.082

6.3 Method: Weibull; p = 100; \(n \in \{250, 500\}\); \(r_{nz} \in \{1\}\)

6.3.1 n = 250

debug_methods(n = 250, p = 100, g = 5, method = "weibull", seed = 40, t_SNR = 4, perf_A = T, nzg = 1) %>%
  system.time()
## **** SETTINGS ****
## n = 250 p = 100 g = 5 method = weibull 
## 
## Generating data...
## [1] "Generating βs..."
## Mean of betas pre-scaling: 0.249922847059995
## Lower range of beta values pre-scaling: -3.52945436691976
## Upper range of beta values pre-scaling: 4.83157475138873
## Method: weibull
## Target SNR: 4
## Epsilon variance: 0.0165
## Mean of betas post-scaling: 2.42861286636753e-19
## Lower range of beta values post-scaling: -0.0112344086238269
## Upper range of beta values post-scaling: 0.0136192146883276
## Achieved SNR: 4

## [1] "Generating X..."
## [1] "Done!"
## [1] "Generating Weibull population lifetable..."
## Range of death: 29.1486009954669
## Range of death: 113.639483758283

## [1] "Done!"

## [1] "Generating data from Weibull lifetable..."
## Nr of valid indices: 2240
## Range of linpred values: -0.2851693 0.2659974 
## [1] "Done!"
## Data generated successfully!
## Dataset size: 500 x 109 
## 
## **** TESTING AFT ****
## AFT model error: Singular hessian; suspicious variable No. TRUE:
##  = 
## Try running with fixed shape
## [1] "Fit failed"
## AFT FAILED!
## 
## **** TESTING BRMS ****
## Start sampling
## Warning in .fun(data = .x1, seed = .x2, init = .x3, adapt_delta = .x4,
## max_treedepth = .x5, : 'stepsize' is deprecated. Please use 'step_size'
## instead.
## Running MCMC with 2 parallel chains, with 2 thread(s) per chain...
## 
## Chain 1 Iteration:    1 / 6000 [  0%]  (Warmup) 
## Chain 2 Iteration:    1 / 6000 [  0%]  (Warmup) 
## Chain 1 Iteration:  100 / 6000 [  1%]  (Warmup) 
## Chain 2 Iteration:  100 / 6000 [  1%]  (Warmup) 
## Chain 2 Iteration:  200 / 6000 [  3%]  (Warmup) 
## Chain 1 Iteration:  200 / 6000 [  3%]  (Warmup) 
## Chain 2 Iteration:  300 / 6000 [  5%]  (Warmup) 
## Chain 1 Iteration:  300 / 6000 [  5%]  (Warmup) 
## Chain 2 Iteration:  400 / 6000 [  6%]  (Warmup) 
## Chain 1 Iteration:  400 / 6000 [  6%]  (Warmup) 
## Chain 2 Iteration:  500 / 6000 [  8%]  (Warmup) 
## Chain 1 Iteration:  500 / 6000 [  8%]  (Warmup) 
## Chain 2 Iteration:  600 / 6000 [ 10%]  (Warmup) 
## Chain 1 Iteration:  600 / 6000 [ 10%]  (Warmup) 
## Chain 2 Iteration:  700 / 6000 [ 11%]  (Warmup) 
## Chain 2 Iteration:  800 / 6000 [ 13%]  (Warmup) 
## Chain 1 Iteration:  700 / 6000 [ 11%]  (Warmup) 
## Chain 1 Iteration:  800 / 6000 [ 13%]  (Warmup) 
## Chain 2 Iteration:  900 / 6000 [ 15%]  (Warmup) 
## Chain 2 Iteration: 1000 / 6000 [ 16%]  (Warmup) 
## Chain 1 Iteration:  900 / 6000 [ 15%]  (Warmup) 
## Chain 2 Iteration: 1100 / 6000 [ 18%]  (Warmup) 
## Chain 1 Iteration: 1000 / 6000 [ 16%]  (Warmup) 
## Chain 2 Iteration: 1200 / 6000 [ 20%]  (Warmup) 
## Chain 1 Iteration: 1100 / 6000 [ 18%]  (Warmup) 
## Chain 2 Iteration: 1300 / 6000 [ 21%]  (Warmup) 
## Chain 2 Iteration: 1400 / 6000 [ 23%]  (Warmup) 
## Chain 1 Iteration: 1200 / 6000 [ 20%]  (Warmup) 
## Chain 1 Iteration: 1300 / 6000 [ 21%]  (Warmup) 
## Chain 2 Iteration: 1500 / 6000 [ 25%]  (Warmup) 
## Chain 1 Iteration: 1400 / 6000 [ 23%]  (Warmup) 
## Chain 2 Iteration: 1600 / 6000 [ 26%]  (Warmup) 
## Chain 1 Iteration: 1500 / 6000 [ 25%]  (Warmup) 
## Chain 2 Iteration: 1700 / 6000 [ 28%]  (Warmup) 
## Chain 1 Iteration: 1600 / 6000 [ 26%]  (Warmup) 
## Chain 2 Iteration: 1800 / 6000 [ 30%]  (Warmup) 
## Chain 1 Iteration: 1700 / 6000 [ 28%]  (Warmup) 
## Chain 2 Iteration: 1900 / 6000 [ 31%]  (Warmup) 
## Chain 1 Iteration: 1800 / 6000 [ 30%]  (Warmup) 
## Chain 2 Iteration: 2000 / 6000 [ 33%]  (Warmup) 
## Chain 1 Iteration: 1900 / 6000 [ 31%]  (Warmup) 
## Chain 1 Iteration: 2000 / 6000 [ 33%]  (Warmup) 
## Chain 2 Iteration: 2100 / 6000 [ 35%]  (Warmup) 
## Chain 1 Iteration: 2100 / 6000 [ 35%]  (Warmup) 
## Chain 2 Iteration: 2200 / 6000 [ 36%]  (Warmup) 
## Chain 1 Iteration: 2200 / 6000 [ 36%]  (Warmup) 
## Chain 2 Iteration: 2300 / 6000 [ 38%]  (Warmup) 
## Chain 2 Iteration: 2400 / 6000 [ 40%]  (Warmup) 
## Chain 1 Iteration: 2300 / 6000 [ 38%]  (Warmup) 
## Chain 2 Iteration: 2500 / 6000 [ 41%]  (Warmup) 
## Chain 1 Iteration: 2400 / 6000 [ 40%]  (Warmup) 
## Chain 2 Iteration: 2600 / 6000 [ 43%]  (Warmup) 
## Chain 1 Iteration: 2500 / 6000 [ 41%]  (Warmup) 
## Chain 2 Iteration: 2700 / 6000 [ 45%]  (Warmup) 
## Chain 1 Iteration: 2600 / 6000 [ 43%]  (Warmup) 
## Chain 2 Iteration: 2800 / 6000 [ 46%]  (Warmup) 
## Chain 1 Iteration: 2700 / 6000 [ 45%]  (Warmup) 
## Chain 2 Iteration: 2900 / 6000 [ 48%]  (Warmup) 
## Chain 1 Iteration: 2800 / 6000 [ 46%]  (Warmup) 
## Chain 2 Iteration: 3000 / 6000 [ 50%]  (Warmup) 
## Chain 2 Iteration: 3001 / 6000 [ 50%]  (Sampling) 
## Chain 1 Iteration: 2900 / 6000 [ 48%]  (Warmup) 
## Chain 1 Iteration: 3000 / 6000 [ 50%]  (Warmup) 
## Chain 2 Iteration: 3100 / 6000 [ 51%]  (Sampling) 
## Chain 1 Iteration: 3001 / 6000 [ 50%]  (Sampling) 
## Chain 1 Iteration: 3100 / 6000 [ 51%]  (Sampling) 
## Chain 2 Iteration: 3200 / 6000 [ 53%]  (Sampling) 
## Chain 1 Iteration: 3200 / 6000 [ 53%]  (Sampling) 
## Chain 1 Iteration: 3300 / 6000 [ 55%]  (Sampling) 
## Chain 2 Iteration: 3300 / 6000 [ 55%]  (Sampling) 
## Chain 1 Iteration: 3400 / 6000 [ 56%]  (Sampling) 
## Chain 1 Iteration: 3500 / 6000 [ 58%]  (Sampling) 
## Chain 2 Iteration: 3400 / 6000 [ 56%]  (Sampling) 
## Chain 1 Iteration: 3600 / 6000 [ 60%]  (Sampling) 
## Chain 1 Iteration: 3700 / 6000 [ 61%]  (Sampling) 
## Chain 2 Iteration: 3500 / 6000 [ 58%]  (Sampling) 
## Chain 1 Iteration: 3800 / 6000 [ 63%]  (Sampling) 
## Chain 1 Iteration: 3900 / 6000 [ 65%]  (Sampling) 
## Chain 2 Iteration: 3600 / 6000 [ 60%]  (Sampling) 
## Chain 1 Iteration: 4000 / 6000 [ 66%]  (Sampling) 
## Chain 2 Iteration: 3700 / 6000 [ 61%]  (Sampling) 
## Chain 1 Iteration: 4100 / 6000 [ 68%]  (Sampling) 
## Chain 1 Iteration: 4200 / 6000 [ 70%]  (Sampling) 
## Chain 2 Iteration: 3800 / 6000 [ 63%]  (Sampling) 
## Chain 1 Iteration: 4300 / 6000 [ 71%]  (Sampling) 
## Chain 1 Iteration: 4400 / 6000 [ 73%]  (Sampling) 
## Chain 2 Iteration: 3900 / 6000 [ 65%]  (Sampling) 
## Chain 1 Iteration: 4500 / 6000 [ 75%]  (Sampling) 
## Chain 1 Iteration: 4600 / 6000 [ 76%]  (Sampling) 
## Chain 2 Iteration: 4000 / 6000 [ 66%]  (Sampling) 
## Chain 1 Iteration: 4700 / 6000 [ 78%]  (Sampling) 
## Chain 1 Iteration: 4800 / 6000 [ 80%]  (Sampling) 
## Chain 2 Iteration: 4100 / 6000 [ 68%]  (Sampling) 
## Chain 1 Iteration: 4900 / 6000 [ 81%]  (Sampling) 
## Chain 2 Iteration: 4200 / 6000 [ 70%]  (Sampling) 
## Chain 1 Iteration: 5000 / 6000 [ 83%]  (Sampling) 
## Chain 1 Iteration: 5100 / 6000 [ 85%]  (Sampling) 
## Chain 2 Iteration: 4300 / 6000 [ 71%]  (Sampling) 
## Chain 1 Iteration: 5200 / 6000 [ 86%]  (Sampling) 
## Chain 1 Iteration: 5300 / 6000 [ 88%]  (Sampling) 
## Chain 2 Iteration: 4400 / 6000 [ 73%]  (Sampling) 
## Chain 1 Iteration: 5400 / 6000 [ 90%]  (Sampling) 
## Chain 1 Iteration: 5500 / 6000 [ 91%]  (Sampling) 
## Chain 2 Iteration: 4500 / 6000 [ 75%]  (Sampling) 
## Chain 1 Iteration: 5600 / 6000 [ 93%]  (Sampling) 
## Chain 1 Iteration: 5700 / 6000 [ 95%]  (Sampling) 
## Chain 2 Iteration: 4600 / 6000 [ 76%]  (Sampling) 
## Chain 1 Iteration: 5800 / 6000 [ 96%]  (Sampling) 
## Chain 1 Iteration: 5900 / 6000 [ 98%]  (Sampling) 
## Chain 2 Iteration: 4700 / 6000 [ 78%]  (Sampling) 
## Chain 1 Iteration: 6000 / 6000 [100%]  (Sampling) 
## Chain 1 finished in 32.8 seconds.
## Chain 2 Iteration: 4800 / 6000 [ 80%]  (Sampling) 
## Chain 2 Iteration: 4900 / 6000 [ 81%]  (Sampling) 
## Chain 2 Iteration: 5000 / 6000 [ 83%]  (Sampling) 
## Chain 2 Iteration: 5100 / 6000 [ 85%]  (Sampling) 
## Chain 2 Iteration: 5200 / 6000 [ 86%]  (Sampling) 
## Chain 2 Iteration: 5300 / 6000 [ 88%]  (Sampling) 
## Chain 2 Iteration: 5400 / 6000 [ 90%]  (Sampling) 
## Chain 2 Iteration: 5500 / 6000 [ 91%]  (Sampling) 
## Chain 2 Iteration: 5600 / 6000 [ 93%]  (Sampling) 
## Chain 2 Iteration: 5700 / 6000 [ 95%]  (Sampling) 
## Chain 2 Iteration: 5800 / 6000 [ 96%]  (Sampling) 
## Chain 2 Iteration: 5900 / 6000 [ 98%]  (Sampling) 
## Chain 2 Iteration: 6000 / 6000 [100%]  (Sampling) 
## Chain 2 finished in 40.1 seconds.
## 
## Both chains finished successfully.
## Mean chain execution time: 36.5 seconds.
## Total execution time: 40.2 seconds.
## Warning: 1 of 6000 (0.0%) transitions ended with a divergence.
## See https://mc-stan.org/misc/warnings for details.
## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.
## Shape estimate: 8.818291 
## Scale estimate 76.27433 
## BRMS SUCCESS!
## BRMS RMSEs: 
## $coef
## [1] 0.007108537
## 
## $mrl
## [1] 6.190739
## 
## $bioage
## [1] 7.873618

## 
## **** COMPARISON ****

## 
## **** SUMMARY COMPARISON ****
## Beta RMSE - AFT: 0 
## Beta RMSE - BRMS: 0.0071 
## MRL RMSE - AFT: 0 
## MRL RMSE - BRMS: 6.1907 
## BioAge RMSE - AFT 0 
## BioAge RMSE - BRMS 7.8736
##    user  system elapsed 
## 128.922  13.561  54.699

6.3.2 n = 500

debug_methods(n = 500, p = 100, g = 5, method = "weibull", seed = 40, t_SNR = 4, perf_A = T, nzg = 1) %>%
  system.time()
## **** SETTINGS ****
## n = 500 p = 100 g = 5 method = weibull 
## 
## Generating data...
## [1] "Generating βs..."
## Mean of betas pre-scaling: 0.249922847059995
## Lower range of beta values pre-scaling: -3.52945436691976
## Upper range of beta values pre-scaling: 4.83157475138873
## Method: weibull
## Target SNR: 4
## Epsilon variance: 0.0165
## Mean of betas post-scaling: 2.42861286636753e-19
## Lower range of beta values post-scaling: -0.0112344086238269
## Upper range of beta values post-scaling: 0.0136192146883276
## Achieved SNR: 4

## [1] "Generating X..."
## [1] "Done!"
## [1] "Generating Weibull population lifetable..."
## Range of death: 29.1486009954669
## Range of death: 113.639483758283

## [1] "Done!"

## [1] "Generating data from Weibull lifetable..."
## Nr of valid indices: 4489
## Range of linpred values: -0.2702238 0.2766213 
## [1] "Done!"
## Data generated successfully!
## Dataset size: 1000 x 109 
## 
## **** TESTING AFT ****
## AFT model error: Singular hessian; suspicious variable No. TRUE:
##  = 
## Try running with fixed shape
## [1] "Fit failed"
## AFT FAILED!
## 
## **** TESTING BRMS ****
## Start sampling
## Warning in .fun(data = .x1, seed = .x2, init = .x3, adapt_delta = .x4,
## max_treedepth = .x5, : 'stepsize' is deprecated. Please use 'step_size'
## instead.
## Running MCMC with 2 parallel chains, with 2 thread(s) per chain...
## 
## Chain 1 Iteration:    1 / 6000 [  0%]  (Warmup) 
## Chain 2 Iteration:    1 / 6000 [  0%]  (Warmup) 
## Chain 1 Iteration:  100 / 6000 [  1%]  (Warmup) 
## Chain 2 Iteration:  100 / 6000 [  1%]  (Warmup) 
## Chain 1 Iteration:  200 / 6000 [  3%]  (Warmup) 
## Chain 1 Iteration:  300 / 6000 [  5%]  (Warmup) 
## Chain 2 Iteration:  200 / 6000 [  3%]  (Warmup) 
## Chain 1 Iteration:  400 / 6000 [  6%]  (Warmup) 
## Chain 2 Iteration:  300 / 6000 [  5%]  (Warmup) 
## Chain 1 Iteration:  500 / 6000 [  8%]  (Warmup) 
## Chain 2 Iteration:  400 / 6000 [  6%]  (Warmup) 
## Chain 1 Iteration:  600 / 6000 [ 10%]  (Warmup) 
## Chain 2 Iteration:  500 / 6000 [  8%]  (Warmup) 
## Chain 1 Iteration:  700 / 6000 [ 11%]  (Warmup) 
## Chain 2 Iteration:  600 / 6000 [ 10%]  (Warmup) 
## Chain 1 Iteration:  800 / 6000 [ 13%]  (Warmup) 
## Chain 2 Iteration:  700 / 6000 [ 11%]  (Warmup) 
## Chain 1 Iteration:  900 / 6000 [ 15%]  (Warmup) 
## Chain 2 Iteration:  800 / 6000 [ 13%]  (Warmup) 
## Chain 1 Iteration: 1000 / 6000 [ 16%]  (Warmup) 
## Chain 2 Iteration:  900 / 6000 [ 15%]  (Warmup) 
## Chain 1 Iteration: 1100 / 6000 [ 18%]  (Warmup) 
## Chain 2 Iteration: 1000 / 6000 [ 16%]  (Warmup) 
## Chain 1 Iteration: 1200 / 6000 [ 20%]  (Warmup) 
## Chain 2 Iteration: 1100 / 6000 [ 18%]  (Warmup) 
## Chain 1 Iteration: 1300 / 6000 [ 21%]  (Warmup) 
## Chain 2 Iteration: 1200 / 6000 [ 20%]  (Warmup) 
## Chain 1 Iteration: 1400 / 6000 [ 23%]  (Warmup) 
## Chain 2 Iteration: 1300 / 6000 [ 21%]  (Warmup) 
## Chain 1 Iteration: 1500 / 6000 [ 25%]  (Warmup) 
## Chain 2 Iteration: 1400 / 6000 [ 23%]  (Warmup) 
## Chain 1 Iteration: 1600 / 6000 [ 26%]  (Warmup) 
## Chain 2 Iteration: 1500 / 6000 [ 25%]  (Warmup) 
## Chain 1 Iteration: 1700 / 6000 [ 28%]  (Warmup) 
## Chain 2 Iteration: 1600 / 6000 [ 26%]  (Warmup) 
## Chain 1 Iteration: 1800 / 6000 [ 30%]  (Warmup) 
## Chain 2 Iteration: 1700 / 6000 [ 28%]  (Warmup) 
## Chain 1 Iteration: 1900 / 6000 [ 31%]  (Warmup) 
## Chain 2 Iteration: 1800 / 6000 [ 30%]  (Warmup) 
## Chain 1 Iteration: 2000 / 6000 [ 33%]  (Warmup) 
## Chain 2 Iteration: 1900 / 6000 [ 31%]  (Warmup) 
## Chain 1 Iteration: 2100 / 6000 [ 35%]  (Warmup) 
## Chain 2 Iteration: 2000 / 6000 [ 33%]  (Warmup) 
## Chain 1 Iteration: 2200 / 6000 [ 36%]  (Warmup) 
## Chain 2 Iteration: 2100 / 6000 [ 35%]  (Warmup) 
## Chain 1 Iteration: 2300 / 6000 [ 38%]  (Warmup) 
## Chain 2 Iteration: 2200 / 6000 [ 36%]  (Warmup) 
## Chain 1 Iteration: 2400 / 6000 [ 40%]  (Warmup) 
## Chain 2 Iteration: 2300 / 6000 [ 38%]  (Warmup) 
## Chain 1 Iteration: 2500 / 6000 [ 41%]  (Warmup) 
## Chain 2 Iteration: 2400 / 6000 [ 40%]  (Warmup) 
## Chain 1 Iteration: 2600 / 6000 [ 43%]  (Warmup) 
## Chain 2 Iteration: 2500 / 6000 [ 41%]  (Warmup) 
## Chain 1 Iteration: 2700 / 6000 [ 45%]  (Warmup) 
## Chain 2 Iteration: 2600 / 6000 [ 43%]  (Warmup) 
## Chain 1 Iteration: 2800 / 6000 [ 46%]  (Warmup) 
## Chain 2 Iteration: 2700 / 6000 [ 45%]  (Warmup) 
## Chain 1 Iteration: 2900 / 6000 [ 48%]  (Warmup) 
## Chain 2 Iteration: 2800 / 6000 [ 46%]  (Warmup) 
## Chain 1 Iteration: 3000 / 6000 [ 50%]  (Warmup) 
## Chain 1 Iteration: 3001 / 6000 [ 50%]  (Sampling) 
## Chain 2 Iteration: 2900 / 6000 [ 48%]  (Warmup) 
## Chain 1 Iteration: 3100 / 6000 [ 51%]  (Sampling) 
## Chain 2 Iteration: 3000 / 6000 [ 50%]  (Warmup) 
## Chain 2 Iteration: 3001 / 6000 [ 50%]  (Sampling) 
## Chain 1 Iteration: 3200 / 6000 [ 53%]  (Sampling) 
## Chain 2 Iteration: 3100 / 6000 [ 51%]  (Sampling) 
## Chain 1 Iteration: 3300 / 6000 [ 55%]  (Sampling) 
## Chain 2 Iteration: 3200 / 6000 [ 53%]  (Sampling) 
## Chain 1 Iteration: 3400 / 6000 [ 56%]  (Sampling) 
## Chain 2 Iteration: 3300 / 6000 [ 55%]  (Sampling) 
## Chain 1 Iteration: 3500 / 6000 [ 58%]  (Sampling) 
## Chain 2 Iteration: 3400 / 6000 [ 56%]  (Sampling) 
## Chain 1 Iteration: 3600 / 6000 [ 60%]  (Sampling) 
## Chain 2 Iteration: 3500 / 6000 [ 58%]  (Sampling) 
## Chain 1 Iteration: 3700 / 6000 [ 61%]  (Sampling) 
## Chain 2 Iteration: 3600 / 6000 [ 60%]  (Sampling) 
## Chain 1 Iteration: 3800 / 6000 [ 63%]  (Sampling) 
## Chain 2 Iteration: 3700 / 6000 [ 61%]  (Sampling) 
## Chain 1 Iteration: 3900 / 6000 [ 65%]  (Sampling) 
## Chain 2 Iteration: 3800 / 6000 [ 63%]  (Sampling) 
## Chain 1 Iteration: 4000 / 6000 [ 66%]  (Sampling) 
## Chain 2 Iteration: 3900 / 6000 [ 65%]  (Sampling) 
## Chain 1 Iteration: 4100 / 6000 [ 68%]  (Sampling) 
## Chain 1 Iteration: 4200 / 6000 [ 70%]  (Sampling) 
## Chain 2 Iteration: 4000 / 6000 [ 66%]  (Sampling) 
## Chain 1 Iteration: 4300 / 6000 [ 71%]  (Sampling) 
## Chain 2 Iteration: 4100 / 6000 [ 68%]  (Sampling) 
## Chain 2 Iteration: 4200 / 6000 [ 70%]  (Sampling) 
## Chain 1 Iteration: 4400 / 6000 [ 73%]  (Sampling) 
## Chain 1 Iteration: 4500 / 6000 [ 75%]  (Sampling) 
## Chain 2 Iteration: 4300 / 6000 [ 71%]  (Sampling) 
## Chain 1 Iteration: 4600 / 6000 [ 76%]  (Sampling) 
## Chain 2 Iteration: 4400 / 6000 [ 73%]  (Sampling) 
## Chain 1 Iteration: 4700 / 6000 [ 78%]  (Sampling) 
## Chain 2 Iteration: 4500 / 6000 [ 75%]  (Sampling) 
## Chain 1 Iteration: 4800 / 6000 [ 80%]  (Sampling) 
## Chain 2 Iteration: 4600 / 6000 [ 76%]  (Sampling) 
## Chain 1 Iteration: 4900 / 6000 [ 81%]  (Sampling) 
## Chain 2 Iteration: 4700 / 6000 [ 78%]  (Sampling) 
## Chain 1 Iteration: 5000 / 6000 [ 83%]  (Sampling) 
## Chain 2 Iteration: 4800 / 6000 [ 80%]  (Sampling) 
## Chain 1 Iteration: 5100 / 6000 [ 85%]  (Sampling) 
## Chain 2 Iteration: 4900 / 6000 [ 81%]  (Sampling) 
## Chain 1 Iteration: 5200 / 6000 [ 86%]  (Sampling) 
## Chain 2 Iteration: 5000 / 6000 [ 83%]  (Sampling) 
## Chain 1 Iteration: 5300 / 6000 [ 88%]  (Sampling) 
## Chain 2 Iteration: 5100 / 6000 [ 85%]  (Sampling) 
## Chain 1 Iteration: 5400 / 6000 [ 90%]  (Sampling) 
## Chain 2 Iteration: 5200 / 6000 [ 86%]  (Sampling) 
## Chain 1 Iteration: 5500 / 6000 [ 91%]  (Sampling) 
## Chain 2 Iteration: 5300 / 6000 [ 88%]  (Sampling) 
## Chain 1 Iteration: 5600 / 6000 [ 93%]  (Sampling) 
## Chain 2 Iteration: 5400 / 6000 [ 90%]  (Sampling) 
## Chain 1 Iteration: 5700 / 6000 [ 95%]  (Sampling) 
## Chain 2 Iteration: 5500 / 6000 [ 91%]  (Sampling) 
## Chain 1 Iteration: 5800 / 6000 [ 96%]  (Sampling) 
## Chain 2 Iteration: 5600 / 6000 [ 93%]  (Sampling) 
## Chain 1 Iteration: 5900 / 6000 [ 98%]  (Sampling) 
## Chain 2 Iteration: 5700 / 6000 [ 95%]  (Sampling) 
## Chain 1 Iteration: 6000 / 6000 [100%]  (Sampling) 
## Chain 1 finished in 92.3 seconds.
## Chain 2 Iteration: 5800 / 6000 [ 96%]  (Sampling) 
## Chain 2 Iteration: 5900 / 6000 [ 98%]  (Sampling) 
## Chain 2 Iteration: 6000 / 6000 [100%]  (Sampling) 
## Chain 2 finished in 94.5 seconds.
## 
## Both chains finished successfully.
## Mean chain execution time: 93.4 seconds.
## Total execution time: 94.6 seconds.
## Warning: 5 of 6000 (0.0%) transitions ended with a divergence.
## See https://mc-stan.org/misc/warnings for details.
## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.
## Shape estimate: 9.720263 
## Scale estimate 76.30104 
## BRMS SUCCESS!
## BRMS RMSEs: 
## $coef
## [1] 0.006483478
## 
## $mrl
## [1] 5.998548
## 
## $bioage
## [1] 7.433909

## 
## **** COMPARISON ****

## 
## **** SUMMARY COMPARISON ****
## Beta RMSE - AFT: 0 
## Beta RMSE - BRMS: 0.0065 
## MRL RMSE - AFT: 0 
## MRL RMSE - BRMS: 5.9985 
## BioAge RMSE - AFT 0 
## BioAge RMSE - BRMS 7.4339
##    user  system elapsed 
## 316.119  27.228 116.807

6.4 Method: Weibull; p = 200; \(n \in \{250, 500\}\); \(r_{nz} \in \{1\}\)

I will not do because takes too long!